Slowly driven sandpile formation with granular mixtures 
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We introduce a one- dimensional sandpile model with iV different particle types and an infinites- 
imal driving rate. The parameters for the model are the A^ 2 critical slopes for one type of particle 
on top of another. The model is trivial when N — 1, but for iV = 2 we observe four broad classes 
of sandpile structure in different regions of the parameter space. We describe and explain the be- 
haviour of each of these classes, giving quantitative analysis wherever possible. The behaviour of 
sandpiles with N > 2 essentially consists of combinations of these four classes. We investigate the 
model's robustness and highlight the key areas that any experiment designed to reproduce these 
results should focus on. 



PACS numbers: 05.40+j, 46.10+z, 64.75+g 



I. INTRODUCTION 

Granular materials display a variety of unusual be- 
haviour not normally associated with either solids or 
liquids pL One such example is the segregation that 
occurs when a mixture of different sized granules is re- 
peatedly shaken, in which the larger particles rise to the 
top Similarly, a granular mixture placed inside a 

rotating cylinder segregates into alternate bands along 
the cylinder's axis ||p[. Segregation in the absence of 
external perturbations has recently been demonstrated 
for a mixture poured between two vertical plates sepa- 
rated by a narrow gap |Q. A sandpile forms in which 
particles of different sizes tend to remain near the top or 
bottom. Moreover, certain mixtures also self-stratify into 
alternating layers parallel to the surface of the pile. It 
is this segregation in sandpiles that this paper seeks to 
address. 

A sandpile is formed by the addition of particles which 
then move over the surface of the pile until finding a 
resting place. Modelling this process is a highly non- 
trivial problem even without the added complication of 
mixtures of particles. There are two basic approaches 
to modelling the surface transport of pure granular me- 
dia. Firstly, Bak et. at [|| introduced a cellular automata 
model as a paradigm of their more general concept of self- 
organised criticality. In this model, sequentially added 
particles can initiate a series of locally-defined topples 
known collectively as an avalanche. The system is said 
to be slowly driven as the timescale for particle in- 
jection is infinitely slower than that of the subsequent 
avalanche H. However, experimental evidence 10-1^] 
disagrees with the model's predicted power-law distri- 
bution of avalanche sizes, exce pt possibly in the limit 
of overdamped particle motion ||13| . A second approach 
treats the sandpile as a continuum with a fluid rolling 
layer interacting with the static bulk of the pile. Analy- 
sis of the field equations of this model appears to give a 
greater correspondence with the experiments flTij . 

Variations of both discrete and continuous approaches 
have been used to try and explain the behaviour observed 



in the sandpile formed by pouring a mixture between 
two plates. A discrete model, in which the particles 
are added in groups and also move downslope in groups, 
exhibits self-stratification but gives no insights into the 
mechanism responsible for the formation of the layers. 
However, extending the continuum rate equations to in- 
corporate mixtures demonstrates the existence of a kink, 
which back-propagates up the slope forming two layers at 
once, giving a possible explanation of the experimental 
results 15 1 . 

In this paper, and in contrast to the experiments JtJ 
and their subsequent analysis |l5|j , we consider a slowly 
driven system in which particles are added individually 
rather than being poured. The entire pile is stable be- 
tween particle additions and there is no formation of a 
rolling layer. The model we have chosen to study is a 
one-dimensional cellular automata with mixed particle 
types, based on the model in js| , which has no avalanch- 
ing in one-dimension. For binary mixtures we observe 
four broad classes of sandpile structure as the relative 
strengths of the interactions between the particles are 
varied. Two of these classes correspond to the self- 
segregation and self-stratification observed in the rolling 
layer case. 

The algorithm for the model is described in Sec.||. The 
four classes of behaviour are described and their evolution 



explained in Sec. [II, and the results for mixtures with 
more than two particle types are also given. In Sec. IV we 



discuss how experiments to observe these classes might 
proceed. 



II. THE MODEL 

The sandpile profile is described by the set of heights 
hi, i > 1, where all the hi are initially set to zero. An in- 
finite wall at i = serves as a lower bound for i, whereas 
there is no upper bound on the values of i. At each time 
step, a particle is chosen from one of N possible types, 
where each particle has the same dimensions and is cho- 
sen with equal probability. A particle of type a £ [l,N] 
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is then added to the top of site 1 and subsequently 
slides to the first site i > 1 that obeys Zi < z a p, where 
Zi = hi — h i+ i and f3 is the type of particle currently on 
the surface at site i. The particle is added to the top 
of site i, hi — ► hi + 1 and another particle can now be 
added to the system. An example of this process is given 
in Fig.| 

The N 2 parameters z a p correspond to the maximum 
slope on which a particle of type a can remain on top of a 
particle of type f3 without sliding off. A pure pile of just 
type a particles has uniform slope z aa , so tan _1 (z QQ ) 
can be identified as the angle of repose. The z a p for 
a ^ (3 are microscopically defined quantities with no ob- 
vious macroscopic counterparts. Generally z a p ^ zp a , 
since the critical slope will depend upon which type of 
particle is moving on top of the pile, for instance types 
a and (3 may have different densities. Due to the com- 
plex nature of granular surface-to-surface interactions it 
is unclear how much of this parameter space corresponds 
to physical reality. The labels for each particle type are 
just dummy variables, so without loss of generality we 
choose to fix zu < Z22 < ■ • • < znn- 
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FIG. 1. Example of sandpile evolution for N — 2, zn = 1, 
212 = 2, Z21 — 3 and 222 — 4. Particles of type 1 are black, 
and those of type 2 are white. For the sandpile shown above, 
a particle of type 1 added to site 1 will eventually come to 
rest on site 4. 



III. RESULTS 

For N = 1 the system reduces to the original model 
which is trivial in one dimension B. For N — 2 and 
zn < Z22 there are four classes of solution, which we 
label I to IV. Each class can be identified according to 
the domain stability of each particle type, which is de- 



fined as follows. A compact region of sites with surface 
composition of particle type a is stable if it can reach a 
uniform slope z aa such that particles of type (3 7^ a will 
slide through the region — that is, if z^ql 2 M . Unsta- 
ble regions can only form in instances where there is no 
incoming flux of particles of the other type, which may 
occur towards the right-hand side of the sandpile. We 
now describe each of the four classes in turn. 

Class I : z 2 i > zn and z 12 > ^22- Neither particle 
type can form stable domains and particles will usually 
come to rest on particles of the other type. We call this 
periodic mixing. An example of such a sandpile is given 
in Fig.^. This behaviour can be demonstrated by the 
following single-site analysis. A site with an uppermost 
particle of type a and slope Zi is represented by (a,Zi). 
Then the transition amplitude for a particle of type (3 
coming to rest on this site, (a, Zj) — » (/?, Zi + 1), is given 
by the particle addition operator P a p{zi) defined by 

P a p{zi) = 0(z af} - z^, (1) 



where 8(x) = 1 for x > and 6(x) — for x < 0. Low 
values for Zi are transitory, but for large Zi added parti- 
cles slide through, possibly coming to rest on site i + 1 
and reducing Zi by 1. Hence the bulk properties of the 
sandpile will be characterised by the action of P in the re- 
gion of large Zj. An example is given in Fig.|| for the case 
zn < Z22 < z\2 < Z21, which clearly shows periodicity for 
Zi near the maximum slope, which in this case is ~ z±2- 
Note that for Z22 > Z21 we do not have strict periodicity 
as type 2 particles are occasionally added consecutively. 




FIG. 2. Example of a class I sandpile with the parameters 
zn — 1, Z12 — 5, 221 = 7 and 222 = 3, which exhibits peri- 
odic mixing. Particles of type 1 are shown as black, particles 
of type 2 are shown as white. The pile given here is small, 
just 35 sites wide at its base, to show the mixing clearly. In 
this and all subsequent cases, simulations have been extended 
up to O(10 5 ) particles without any observed deviation from 
the characteristic behaviour. Note that we have rescaled the 
t/-axis to give a roughly square picture. 
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FIG. 3. Phase portrait of the particle addition operator 
Paffizi) in the case zn < 222 < 212 < Z21, where each ar- 
row corresponds to P a g(zi) = 1, that is, an allowed particle 
addition, a denotes the type of particle on top and z% is the 
slope. 

Class II : z 2 i > zn and z 12 < 222- A stable domain of 
type 2 particles with uniform slope Z22 builds up to the 
left of the sandpile, and type 1 particles slide through to 
the right forming a domain with uniform slope Z\x- The 
result is discrete self-segregation, as the example in Fig.|] 
demonstrates. The boundary between the two domains 
moves to the right as particles of type 2 come to rest on 
sites with slope zn, so if 212 > zn some periodic mixing 
may occur. If 221 > 212, the mixing remains confined to a 
narrow layer at the boundary, but if 221 < 212 it expands 
to the bulk of the pile, giving a sandpile similar in ap- 
pearance to case of the periodic mixing described previ- 
ously. The self-segregation and boundary behaviour can 
also be seen from the phase portraits of the particle ad- 
dition operator P a fi(zi). We justify referring to the case 
211 < 221 < 212 < 222 as discrete self-segregation by the 
observation that the region of periodic mixing collapses 
to a much narrower layer under most small alterations in 
the dynamical rules. We will return to this point later 
when we discuss the robustness of the model. 

It is possible to construct a global solution of this class 
of sandpiles. Define Lb and L to be position of the 
boundary and the right-hand edge of the sandpile re- 
spectively, so that fiL > 0, fiL+i = 0, all sites 1 < i < Lb 
have type 2 particle on top and all sites Lb < i < L have 
type 1 particles on top. A type 2 particle added to i = Lb 
will reduce the slope of site Lb — 1 by 1, so the next type 2 
particle added will stop at Lb — 1, then Lb — 2, Lb — 3, 
and so on. Similarly, type 1 particles will be added to L, 
L — 1, L — 2, . . . , Lb + 1, in that order. L will move to 
the right by one step when increases from to 21 1, so 
in the continuum limit 
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L- L, 
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where the time scale has been normalised to one par- 
ticle addition of either type per unit time. The slope 



at Lb increases with each type 2 particle that stops on 
Lb and decreases with each type 1 particle that stops 
on Lb + 1, and since Lb increases by one whenever the 
slope increases from 21 1 to 222 we get 
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dt 
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(3) 



A steadily evolving sandpile corresponds to a constant 
ratio Lb/L, so 
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dt 
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Using this together with (||) and (||), we find 
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222/211 - 1 

giving the average slope of the entire pile as 



L, 



L, 



1 - — J 211 + —^22 = V Z H Z 22, 



in agreement with simulations |16| 
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FIG. 4. Example of a class II sandpile with zii = 1, 
212 = 3, Z21 = 5 and 222 = 7, demonstrating discrete 
self-segregation. Particles of type 1 are black, those of type 2 
are white. The base of the pile is 100 sites wide. 

Class III : 221 < 2n and 212 > 222- It might be ex- 
pected that self-segregation will also occur here, this time 
with type 1 particles to the left of the boundary. How- 
ever, since z\\ < 222 the type 2 particles to the right of 
the boundary form steeper slopes than the type 1 par- 
ticles to the left, and so now the boundary moves ups- 
lope rather than downslope, creating a double layer of 
2's on top of l's. Once the boundary reaches the left- 
hand wall, a thin layer of l's quickly cover the surface of 
the whole sandpile, the boundary returns to the bottom 
and starts propagating upwards once more. This mov- 
ing interface corresponds to the kink described in the 
rolling layer case pljLq], and, indeed, the region of pa- 
rameter space in which it occurs is the same. However, 
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there are 3 significant differences between the nature of 
the self-stratification formed by this slowly-driven pro- 
cess and that formed by rolling layers, (i) The slope of 
the layers varies between z\\ and Z22, as opposed to the 
uniform slope of 222 observed in the rolling layer case, 
(ii) The layers are narrower, typically just one particle 
wide, (iii) The rate at which the interface moves is no 
longer a constant but varies according to statistical fluc- 
tuations in the types of incoming particles. It is even pos- 
sible for the interface to stop moving altogether, resulting 
in a vertical build-up of particles that is reminiscent of 
a miniature self-segregated sandpile. However, this state 
is unstable and the interface will eventually start moving 
again, either upslope to continue the layering process or 
quickly downslope to the bottom of the pile. Thus, lay- 
ers can now start and stop in the bulk of the pile. An 
example of a self-stratified sandpile is given in Fig.[|. 




FIG. 5. Example of a class III sandpile with zu = 5, 
Z12 = 15, 221 = 1 and 222 — 10, demonstrating 
self-stratification. Type 1 particles are black and type 2 par- 
ticles are white. The base of the pile is 100 sites wide. 

Class IV : Z21 < Z\\ and zvi < 222- In this final class, 
alternating stable domains of l's and 2's form parallel 
vertical bands. An example of this vertical stratifica- 
tion is given in Fig.^j. The phase portrait of P a /3(zi) 
for this class demonstrates the separation of the parti- 
cle types, but further analysis requires some knowledge 
of the global solution. Once the slope of the site at the 
right-hand edge of any domain decreases by one, a single 
layer of particles of similar type back-propagates to the 
left-hand edge, when the adjacent domain will undergo a 
similar process. Thus the bulk of the sandpile builds up a 
layer at a time in this piecewise fashion, from the bottom 
of the pile to the top. The process that initially generates 
each domain depends upon the interaction between the 
bulk of the sandpile and the qualitatively different end 
region at the far right-hand side. In the end region, a 
layer of 2's of thickness w z\\ — Z21 + 1 back-propagates 
to the first domain of type 2 particles. The domain then 
broadens to the right, and the process is repeated. Alter- 
natively, if 212 > Z21 the back-propagating layer of 2's can 
be stopped prematurely by an incoming flux of l's, re- 



sulting in the formation of a new stable domain of type 1 
particles. Since the particles that flow into the end re- 
gion are precisely those that did not stop in the bulk of 
the pile, the number of domains of each different type 
tends to remain equal, although the domains themselves 
get broader as the end region expands. If z%2 < 221 the 
formation of new type 1 domains can still occur when the 
sandpile is small is small due to statistical fluctuations, 
but this rarely occurs for larger piles which take on a 
self-segregated appearance. 




FIG. 6. Example of a class IV sandpile with 211 = 5, 
212 = 3, 221 = 1 and 222 = 7, demonstrating vertical stratifi- 
cation. Type 1 particles are white, type 2 particles are white 
and the base of the pile is 100 sites wide. 



A diagram of the parameter space showing the regions 
in which each of these classes occurs is given in Fig.0. The 
borders between these regions correspond to Z21 = zii 
or Z12 — 222, when domain stability is not well defined. 
In these cases the sandpile behaviour is either indeter- 
minate between the two classes in question or just re- 
duces to random mixing. For z\\ = 222 periodic mixing 
and vertical stratification are unaffected but there is no 
longer any distinction between discrete self-segregation 
and self-stratification. Instead, these two classes are re- 
placed with a hybrid class that exhibits self-segregation 
with a broad, layered boundary. 



II. Discrete 
Self-segregation 



IV. Vertical 
Stratification 



I. Periodic 
Mixing 



III. Self-stratification 



FIG. 7. Schematic diagram of 2 2 i — 2n versus 212 — 222 for 
N = 2 and 211 < 222, showing where each class of sandpile so- 
lution applies. 
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Each class is said to be robust if its existence is insensi- 
tive to the exact choice of dynamical rules. It is possible 
to vary the volumes of each particle type added or to 
introduce an open or closed right-hand boundary condi- 
tion without any significant alteration in the resultant 
sandpile. Similarly, initially adding the particles over a 
range of sites does not affect the sandpile to the right 
of the range, and introducing annealed disorder to the 
Zap just increases the noise. More significant is the effect 
of averaging the z a p over adjacent pairs of sites, corre- 
sponding perhaps to the nestling of the upper particle 
in between the two lower ones, which destroys vertical 
stratification and instead gives random mixing or discrete 
self-segregation. Crudely modelling inertia by allowing 
the moving particle to stochastically drift a short dis- 
tance further than normal distorts vertical stratification 
and replaces layering and periodic mixing with random 
mixing. We conclude that the model is robust except 
when we include displacements in the particle's horizon- 
tal motion. 

For N > 2 the parameter space for z a p becomes too 
large to explore systematically. The situation improves 
somewhat if only domain stability is considered, but this 
still leaves 2 N ( N ~ 1 ' > possible combinations, so we have 
limited ourselves to a brief survey of all these cases for 
N = 3 and a representative sample for N = 4. The re- 
sultant sandpiles are essentially just combinations of the 
four classes identified for N = 2, the only significant new 
feature being that periodic mixing can now occur with 
periodicity < N. Any particle separation and the allowed 
orders of periodic mixing can be predicted in each case 
by extending the particle addition operator P a p{zi) to 
include a, (3 <E [1, JV] . As before, there is no generic way 
of constructing global solutions but in many cases known 
solutions for N — 1 particle systems can be used instead. 
This is possible when particle types a and (5 would by 
themselves periodically mix, when to good approxima- 
tion they can be replaced by a single particle type a' 
which alternates between the two. A domain is stable 
to a' only if it is stable to both a and /3, a domain of 
type a 1 particles is stable to another type only if both a 
and /3 are, and z a i a i = min(z Q ^, zp a ). This reduced sys- 
tem usually exhibits the correct qualitative structure of 
the original but significantly underestimates the amount 
of noise. 



IV. DISCUSSION 

The discrete model studied here is perhaps the simplest 
conceivable model describing slowly driven sandpile for- 
mation with granular mixtures. Nonetheless it exhibits 
a wide variety of non-trivial behaviour in one-dimension, 
and we can only suppose that it will continue to do so 
in higher dimensions. The behaviour of sandpiles for bi- 
nary mixtures falls into one of four classes, two of which 
have known counterparts in the rolling layer model |@,|f5| . 



They are also much more susceptible to statistical fluc- 
tuations in the order of particle types added. 

We have recently initiated a series of experiments in 
an attempt to reproduce these classes with real granular 
materials JT^ ]. On first inspection, there may appear to 
be little hope that such a simple model could describe 
real sandpile mixtures. For instance, the possibility of 
particles bouncing or dislodging surface material has not 
been catered for. However, our numerical investigations 
into the robustness of the model leads us to suppose that 
agreement might be possible in the limit of overdamped 
particle motion. It is also important to realise the limit 
of infinitely slow driving as closely as possible, to min- 
imise the probability of a rolling layer forming, because 
even a thin rolling layer would displace surface material 
along the length of the slope. This would interfere with 
the formation of both vertical stratification and periodic 
mixing. 

In summary, sandpile formation by granular mixtures 
exhibits a greater diversity of behaviour in the slowly 
driven limit than in the rolling layer case, at least numeri- 
cally. If toppling were included into this model |T^] , more 
than one particle would be able to move simultaneously 
in the form of an avalanche. Similarly, allowing mov- 
ing particles to dislodge surface material in some manner 
might allow for something akin to a rolling layer to form. 
It would be interesting to see if the system diversity was 
reduced in either of these two cases. 



ACKNOWLEDGEMENT 

We would like to thank Peter Hobson for useful discus- 
sions. 



Electronic address: David.Head@brunel.ac.uk 
Electronic address: G.J.Rodgers@brunel.ac.uk 

[1] H.M.Jaeger and S.R.Nagel, Science 255, 1523(1992). 

[2] A.Rosato, K.J.Strandburg, F.Prinz and R.H.Swendsen, 
Phys.Rev.Lett. 58, 1038(1987). 

[3] J. B. Knight, H.M.Jaeger and S.R.Nagel, Phys.Rev.Lett. 
70, 3728(1993). 

[4] W.Cooke, S.Warr, J.M.Huntley and R.C.Bali, 
Pbys.Rev.E 53, 2812(1996). 

[5] K.M.Hill and J.Kakalios, Phys.Rev.E 49, R3610(1994). 

[6] O.Zik, D.Levine, S.G.Lipson, S.Shtrikman and 
J.Stavans, Phys.Rev.Lett. 73, 644(1994). 

[7] H.A.Makse, S.Havlin, P.R.King and H.E.Stanley, sub- 
mitted to Nature. 

[8] P.Bak, C.Tang and K.Wiesenfeld, Phys.Rev.A 38, 
364(1988). 

[9] G.Grinstein in Scale Invariance, Interfaces and Non- 
Equilibrium Dynamics, Vol. 344 of NATO Advanced 



5 



Study Institute, Series B:Physics, edited by A.McKane 
et.al. (Plenum, New York, 1995). 
[10] H.M.Jaeger, C.Liu and S.R.Nagel, Phys.Rev.Lett. 62, 
40(1989). 

[11] G. A. Held, D.H.Solina, 
II, D.T.Keane, W.J.Haag, P.M.Horn and G.Grinstein, 
Phys.Rev.Lett. 65, 1120(1990). 

[12] P.Evesque and J.Rajchenbach, Phys.Rev.Lett. 62, 
44(1989). 

[13] V.Frette, K.Christensen, A.Malthe-S0renssen, J.Feder, 
T.j0ssang and P.Meakin, Nature 379, 49(1996). 



[14] J.-P.Bouchad, M.E.Cates, J.R.Prakash and 
S.F.Edwards, Phys.Rev.Lett. 74, 1982(1995). 

[15] H.A.Makse, P.Cizeau and H.E.Stanley, submitted to 
PRE or PRL. 

[16] This result can also be derived by a geometric argument 
if it is assumed that both domains sweep out equal ar- 
eas in time, corresponding to the equal rates of particle 
addition. 

[17] D. A. Head, P.R.Hobson and G.J.Rodgers, in preparation. 
[18] L.P.Kadanoff, S.R.Nagel, L.Wu and S.Zhou, Phys.Rev.A 
39 6524(1989). 



G 



